Cite this Article
Gu Yun-Feng, Zhu Liu-Tong, Wu Xiao-Li. Polarization resolved analysis of phonon transport in a multi-terminal system. Chinese Physics B, 2019, 28(12): 124401  
Permissions
Polarization resolved analysis of phonon transport in a multi-terminal system
Gu Yun-Feng †, Zhu Liu-Tong, Wu Xiao-Li
College of Electronic and Mechanical Engineering, Nanjing Forestry University, Nanjing, 210037, China

 

† Corresponding author. E-mail: gu_yunfeng@sina.com

Project supported by the National Natural Science Foundation of China (Grant No. 51376094) and Jiangsu Overseas Visiting Scholar Program for University Prominent Young & Middle-aged Teachers and Presidents, China.

Abstract

The atomistic Green’s function method is improved to compute the polarization resolved phonon transport in a multi-terminal system. Based on the recent developments in literature, the algorithm is simplified. The complex phonon band structure of a semi-infinite periodic terminal is obtained by the generalized eigenvalue equation. Then both the surface Green’s function and phonon group velocity in the terminal are determined from the wave modes propagating away from the scattering region along the terminal. With these key ingredients, the individual phonon mode transmittance between the terminals can be calculated. The feasibility and validity of the method are demonstrated by the chain example compared with the wave packet method, and an example of graphene nanojunction with three terminals.

PACS: 44.10.+i;63.20.D-;63.22.-m;81.05.ue
Keyword:Green’s function method;phonon transport;multi-terminal systems
1. Introduction

The atomistic Green’s function (AGF) method is a powerful tool to explore the phonon transport mechanism in the semiconductor and dielectric nanostructures.[1] Since Mingo and Yang[2] developed this method, extended from the nonequilibrium Green’s function (NEGF) method used to study the electron transport,[3] many important improvements have been made. The AGF method is good at analyzing phonon transport across interfaces. In the traditional AGF method, the calculated phonon transmission function is a summation of the transmittances of all the phonon modes at a given frequency. However, it is very helpful to understand and engineer the thermal transport mechanism on the level of individual phonon polarizations.[4–12] Huang et al.[7] computed the polarization resolved transmittance by using a decomposed self-energy calculated from the density of states matrix for each polarization, but with obscure dependence on the phonon dispersion. Ong et al.[8] extended the AGF method to calculate the transmittance of individual phonon modes across the graphene–hexagonal boron nitride interface based on the concept of the Bloch matrix. The Bloch matrix can be found from the surface Green’s function of the materials on either side of the interface. The phonon-eigenspectrum-based AGF method was developed by Sadasivam et al.[10] for studying the transmission and conversion of individual phonon modes at the interface. For now, the polarization resolved AGF method has been mainly used for the two-terminal systems, for example, the interface between two different materials, and the algorithm is complicated. Many researches show that the multi-terminal systems have some interesting properties, such as phonon Hall effect,[13] thermal rectification,[14] and thermoelectricity,[15] to name a few.

In this paper, we report a multi-terminal polarization resolved AGF method. The scheme is relatively simple to apply, capable of exploring individual phonon mode transmittance between the terminals through a scattering region with arbitrary geometry. The paper is organized as follows. In Section 2, the method used in the calculation is given, including a brief introduction of the traditional AGF method for multi-terminal systems, the surface Green’s function calculation based on the generalized eigenvalue equation, and the polarization resolved transmittance. In Section 3, we give two numerical examples. The first chain example is employed to illustrate the feasibility and validity of the presented method, along with the wave packet (WP) simulation results as a comparison. The second example is about a graphene nanojunction with three terminals. Finally, some conclusions are drawn from the present study.

2. Methods

The model configuration discussed in this paper is a multi-terminal system. Figure 1 shows an example with three terminals. The central part is a scattering region, which is attached with several semi-infinite terminals, serving as the heat reservoirs. The objective of this paper is to explore the way to study the phonon ballistic transport across the scattering region between the terminals.

Fig. 1. An example of a multi-terminal system. The numbers in the terminals denote the respective unit layer indices.

The traditional AGF method can give the total phonon transmittance between terminals. The improved one proposed herein can furthermore study the transmittance between modes in different terminals. As discussed below, both traditional and improved AGF methods can be based on a same generalized eigenvalue equation associated with the periodic properties of the terminal.

2.1. Total transmittance

The dynamical equation for the system, as the starting point, can be expressed as[16]

where I is the identity matrix, and K is the dynamical matrix. ω is the angular frequency. u is the mass-normalized displacement defined as the atomic displacement vector multiplied by the square root of the corresponding atomic mass.[5] For a multi-terminal system, as an example shown in Fig. 1, the dynamical matrix in Eq. (1) can further be divided into[17]
where Km (m = 1,2,…,N) and KS stand for the dynamical matrices of the N terminals and the scattering region (S), respectively. The coupling between the terminals and the scattering region is indicated by in Eq. (2), where the superscript “T” herein represents the matrix transpose. Because the terminals are semi-infinite, the size of Km or Km,S is infinite.

According to the AGF method, the Green’s function for the scattering region has the form[17]

where , and 0+ is the broadening constant, which is a small positive number standing for the phonon energy dissipation in the terminals.[18] i and I are respectively the unit of imaginary number and the identity matrix. In Eq. (3),
where Σm is the self-energy matrix of the m-th terminal, serving as the boundary conditions of terminal applied to the scattering region boundaries.[3] gm is the surface Green’s functions of the terminals. The terminals have periodic layer structures, and the unit layers are numbered in sequence according to the distance from the scattering region as shown in Fig. 1. Because the scattering region atoms have interactions with only the atoms of the 0-th layer in terminals, has a finite number of nonzero elements. As a result, only the submatrix in gm corresponding to atoms in the 0-{th} layer is needed in the calculation of Σm by using Eq. (4). Furthermore, the self-energy matrix Σm has the same size of the dynamical matrix of scattering region KS, and has only a small number of nonzero elements related to the 0-{th} layer atoms in the m-th terminal.[19] It should be noted that there is no terminal–terminal coupling in this paper, and different terminals have interaction with different parts of the scattering region as shown in Fig. 1. This can be achieved by taking a large enough scattering region. Therefore, the location of the nonzero elements in the self-energy matrix Σm is different for different terminals. The simple summation of self-energy matrices Σ in Eq. (4) stands for the effect due to all the terminals on the scattering region. Then we can project the dynamics of the infinite system shown in Fig. 1 to the finite scattering region via the self-energies.[20] Inserting Σ in Eq. (3) results in the Green’s function for the scattering region GS.

The self-energy matrix Σm also can be used to construct the broadening function written as

Based on GS and Γm, the traditional AGF method gives the total transmittance between the terminals[18]
The superscript “*” represents the complex conjugate transpose operation. It should be noted that the transmission function in Eq. (6) is the summation of the transmittances of all modes. The purpose of this paper is to figure out the contribution to the total transmittance τmn from each mode.

2.2. Surface Green’s function based on the generalized eigenvalue equation

The key ingredient in the AGF method is to calculate the surface Green’s function . As described above, the surface Green’s function is a crucial step to get self-energy matrix Σm in Eq. (4). Various methods have been used, and the most accurate and fastest one is based on the solution of the generalized eigenvalue equation.[5] This algorithm takes into account of the translational symmetry of the semi-infinte terminal with periodic layer structures, and can give the phonon dispersion relations in the terminal.

As a submatrix of K in Eq. (2), the dynamical matrix of the m-th terminal with identical unit layers numbered in sequence as shown in Fig. 1 has the form

where is the onsite dynamical matrix of unit layer in the m-th terminal, and denotes the connection between unit layers. The unit layer must be large enough so that only the neighbor unit layers have interactions. Based on the dynamical matrix of the terminal, the phonon dispersion of the terminal can be described by the generalized eigenvalue equation[5,10]
where denotes the phonon eigenvector corresponding to the j-th unit layer in the m-th terminal. Because of the periodicity, the eigenvectors satisfy Bloch theorem, so that defined in the second line of Eq. (8), where the Bloch phase factor λm = ei qam. In this factor, q is the wave number, and am denotes the period of the unit layer along the terminal direction.

There are two groups of solutions for Eq. (8), which indicate waves propagating in opposite directions in the terminal. There is a convenient way to determine which group a wave mode belongs to. By replacing ω with , where 0+ is a small positive perturbation, the outgoing waves propagating away from the scattering region along the terminal will have . All the normalized vectors of the outgoing waves are combined together to form a matrix , and the corresponding eigenvalue matrix is a diagonal matrix , where p is the number of waves at frequency ω. It should be noted that there are two kinds of waves , i.e., propagating modes and evanescent modes. The propagating modes always have a real wave number and propagate stably in the terminal. The evanescent modes, however, will exponentially decay as going away from the scattering region because of the complex wave number , and its imaginary part is positive.

Based on the outgoing wave solutions of the generalized eigenvalue equation with , the surface Green’s function of the m-th terminal is given by[5]

where P is the propagator[6] used to describe the propagation of the phonon waves along the unit layers in the terminal. In the derivation of Eq. (9), is used. In order to make sure that is finite in the limit s → ∞, only the outgoing wave modes are allowed in the propagator P.[5]

2.3. Polarization resolved transmittance

The solutions of the generalized eigenvalue equation (8) not only can be used to get the surface Green’s function of terminal in Eq. (9), but also is helpful to improve the AGF method to analyze the transmittance in terms of phonon polarization, frequency, and wave vector. The remaining piece of ingredient needed in the improved AGF method is the group velocity matrix of the m-th terminal, which is defined as[8]

It should be noted that no contribution to heat flux comes from the evanescent modes which have zero group velocity calculated from Eq. (10).

Now we can proceed to calculate the transmittance between individual phonon channels in different terminals based on the matrix[8]

Let represent the (k,l)-th element of matrix tmn. Then the transmittance between the l-th phonon mode in the n-th terminal and the k-th phonon mode in the m-th terminal is defined as the square modulus
The proposed method herein is straightforward and has a clear physical picture. The process of phonon propagation at a given frequency can be described by the dynamical equation (1) in the frequency space. Equation (1) gives a collective vibration picture. Because the terminal serves as a wave guide with a periodic layer structure, the collective vibration in the terminal can be determined by Eq. (8), as the requirement of the Bloch theorem. The corresponding wave periodicity in space is characterized by the wave vector along the terminal.[21] For the terminal as a one-dimensional system, it is simple to use the scalar quantity, i.e., the wave number. At a given frequency, there are many propagating modes identified by the real wave numbers in the terminals. For any terminal, once , , and are obtained, we can follow a uniform and standard procedure to get all the allowed wave modes by using Eq. (8). The transmittance between these modes in different terminals is determined by Eq. (12).

3. Numerical examples
3.1. Examples of linear atomic chain

The toy examples based on the atomic chain with one degree of freedom per atom, as illustrated in Fig. 2, are employed to illustrate the improved AGF method introduced above. These examples are simple enough for readers to reproduce the results. The central part of the chains is the scattering region having two atoms with masses m and 2m, respectively. In chain (a) and chain (b), there are two semi-infinite chains attached to the scattering region as terminals. The terminal T1 has a one-atom unit cell with the atomic mass of m. Both the terminals T2 and T3 are diatomic chains, and the two atoms in the unit cell have masses 1.2m and 0.8m (1.2m and 1.8m) for the terminal T2 (T3), respectively. The chain (c) in Fig. 2 has three terminals. The inter-atomic spring constant between adjacent atoms is g. All atoms are separated by d when at rest. The period of the unit layer in the terminal, therefore, is d and 2d for T1 and T2 (T3), respectively.

Fig. 2. Schematic of a junction with semi-infinite linear atomic chains as terminals. The scattering region comprises of two atoms with masses m and 2m, respectively. The chain (a) and chain (b) have two terminals, and the chain (c) has three terminals. There are three kinds of terminals: T1, T2, and T3, denoted by red, blue, and green lines, respectively. The numbers in the terminals denote the respective unit layer indices.

The three semi-infinite chain terminals T1, T2, and T3 in Fig. 2 have simple phonon dispersion relations illustrated in Fig. 3. The lines stand for the analytical results,[22] and the red, blue, and green lines correspond to T1, T2, and T3, respectively. The dots represent the outgoing wave modes calculated by using the generalized eigenvalue equation (8) with . As shown in Fig. 3, half modes in the Brillouin zone belong to the outgoing wave modes which have positive group velocities. The diatomic chain T2 or T3 has a frequency gap between the bottom of the optical branch and the top of the acoustic branch. Phonon modes in the frequency gap are evanescent modes which have imaginary wave numbers.[23] These evanescent modes, not shown in Fig. 3, can not propagate along the chain stably.

Fig. 3. Phonon dispersion relations for the terminal chains T1, T2, and T3 shown in Fig. 2. The lines and dots stand for the analytical results[22] and the calculated outgoing waves, respectively. is the characteristic frequency scale.

When a specific phonon mode in one terminal is incident on the scattering region, a part of its energy is transferred to the phonon modes in other terminals. The improved AGF can be used to figure out the transmission coefficients by using Eq. (12). The results are illustrated by the lines in Fig. 4 for the three chain structures discussed above. In fact, the physical meaning of the transmittance can be revealed more intuitively by a completely different method, i.e., the WP method. As a comparison, and in order to give a solid foundation for the AGF calculation, the WP results are also given as circle or triangle symbols in Fig. 4.

Fig. 4. Frequency dependence of phonon transmittance for (a) the two-terminal chains and (b) the three-terminal chain shown in Fig. 2. The lines and symbols stand for the improved AGF results and WP results, respectively.

In the WP method, a WP is generated in T1 as the initial conditions. The WP can be regarded as a phonon. Constructed on the basis of the phonon dispersion relation denoted by the red line in Fig. 3, the WP can be expressed as[24]

Here, u(xn,t;xc,qd) is the time dependent displacement of the atom at location xn = nd in T1. The WP is initially centered at position xc, and has dominant wave vector qd. The amplitude is U = 0.01d. The spatial width of the WP is set to be W = 16d. The group velocity is vg = d ω/dq, calculated based on the phonon dispersion relation of T1 in Fig. 3. The initial atom velocity is defined as v(xn,t;xc,qd) = d u(xn,t;xc,qd)/dt. The WP once is launched in T1, a part of its energy is transferred into another terminal as a new WP. Then the transmission coefficient in the WP method is defined as
where E and Etr are the total energies (including the potential and kinetic energies) of the incident and transmitted WP, respectively. The WP results are denoted by circle and triangle symbols in Fig. 4, and agree well with the AGF results.

The chain (a) in Fig. 2 is the same one studied by Ong,[11] and the same result, denoted by the solid blue line in Fig. 4(a), is obtained. When ω < 1.29ω0, the transmittance between the acoustic phonons in T1 and T2 decreases slowly as the frequency rises, and at last suddenly drops to zero because of the phonon band gap in T2 as shown in Fig. 3. There is a mode conversion between the acoustic phonons in T1 and the optical phonon in T2 with frequency in the range from about 1.58ω0 to 2.04ω0. The chain (b) and chain (a) in Fig. 2 have a similar structure, but different atomic mass in the right terminal. The transmission coefficients of the two chains, therefore, have similar frequency dependence, denoted by the lines in Fig. 4(a).

The chain (c) in Fig. 2 has three terminals. τ12, τ13, and τ23 in Fig. 4(b) represent the phonon transmission coefficients between T1 and T2, T1 and T3, T2 and T3, respectively. In the WP simulation, after scattering with the scattering region, the WP in T1 is clearly observed to be split into two WPs: one is in T2 and the other is in T3. Therefore, it is reasonable that in the low frequency range, say ω < 0.8 ω0, τ12 shown in Fig. 4(b) is significantly lower than the chain (a) result in Fig. 4(a). The same phenomenon is found between τ13 in Fig. 4(b) and the chain (b) result in Fig. 4(a).

As the frequency further rises, τ13 drops rapidly to zero because of the phonon band gap in T3, while τ12 increases significantly. The maximum value of τ12 is close to unity, much higher than the transmittance in chain (a) shown in Fig. 4(a). Furthermore, when the frequency is above the maximum allowed frequency in T3, the larger τ12 reveals that the mode conversion between the acoustic phonons in T1 and the optical phonon in T2 is enhanced. It is interesting that the existence of T3 is helpful for phonon energy exchange between T1 and T2 when the frequency is in the phonon band gap of T3. This extra-large τ13 is also observed when the frequency is in the phonon band gap of T2. The simple theoretical three-terminal model illustrated in Fig. 2(c) is thus capable of guiding phonons into different terminals according to the frequency.

3.2. Examples of graphene nanoribbon

A more complicated example is illustrated in Fig. 5, which is a graphene nanojunction with three same zigzag graphene nanoribbons as terminals. Nz denotes the number of zigzag lines, and Nz = 12 in this example. The fourth nearest-neighbor approximation[25] is used to build the dynamical matrix. In order to avoid the interaction beyond neighboring unit layers in the terminal, the period of the unit layer a has a double size of the primitive cell in the zigzag graphene nanoribbon.

Fig. 5. Schematic diagram of a graphene nanojunction with three terminals labeled as T1, T2, and T3, respectively. Only one unit layer indicated by red dots is shown for each terminal. The period of the unit layer along the terminal direction is a. The width of the zigzag graphene nanoribbon is denoted by Nz which is the number of zigzag lines.

The green or grey lines illustrated in Fig. 6 stand for the phonon branches of the zigzag graphene nanoribbon terminal in Fig. 5. Because of the large period of the unit layer, the size of the Brillouin zone is 2π/a in this simulation due to the phonon folding. The dots represent the outgoing wave modes calculated by using the generalized eigenvalue equation (8) with . There are four acoustic branches in the graphene nanoribbon:[26] the longitudinal (LA) branch, the transverse (TA) branch, the flexural (ZA) branch, and the fourth acoustic (ZA4) branch having a rotational symmetry around a central axis along the nanoribbon, which are represented by the blue, red, black, and magenta dots in Fig. 6, respectively.

Fig. 6. Phonon dispersion of zigzag graphene nanoribbon with Nz = 12. The size of the Brillouin zone is 2π/a. The dots stand for the calculated outgoing waves. The blue, red, black, and magenta dots represent the LA, TA, ZA, and ZA4 branches, respectively.

The total transmittance between the terminals of the graphene nanojunction as a function of phonon frequency is illustrated in Fig. 7(a). As a comparison, the total transmittance τpristine of the pristine zigzag graphene nanoribbon is given as a dash line. τpristine takes on integer values depending on the number of phonon modes, and the value is 4 at the lowest frequency due to the four acoustic branches. τpristine gives the upper limit of the total transmittance, because all the modes are taken into account without scattering.

Fig. 7. (a) The total transmittance between terminals, (b) the polarized transmittance τ12 and (c) τ13 of four acoustic phonon branches in T1. The dash line in (a) is for the pristine graphene nanoribbon as a comparison.

The blue, red, and green lines respectively represent the total transmittance between terminals τ12, τ13, and τ23, which are calculated by using Eq. (6) of the traditional AGF method. The blue, red, and green dots are obtained by summation of mode contribution by using Eq. (12) of the improved AGF method. The dots match well with the lines. As shown in Fig. 5, T2 can be regarded as a perturbation to the straight graphene nanoribbon. Therefore, τ13 has the highest value. The difference between τ12 and τ23 is the result of the asymmetric structure.

The improved AGF method is capable of revealing the polarized transmittance based on Eq. (12). For a specific phonon mode in one terminal, a part of its energy will transfer and turn into various phonon modes in other terminals. The transmittance of the four acoustic branches in T1 is given in Figs. 7(b) and 7(c). It appears that the acoustic phonons in T1 prefer to propagate along the straight nanoribbon into T3 rather than into T2. As the phonon frequency rises, for the ZA modes, τ13 increases significantly and approaches to unity, while τ12 reduces and finally tends to zero.

4. Conclusion

An improved AGF method is developed for revealing the polarized phonon transmittance between terminals in the multi-terminal system, and the algorithm complexity is dramatically decreased. The system has a scattering region with arbitrary geometry, attached with several terminals having a semi-infinite periodic structure. The phonon complex band structure in the terminal is obtained by the generalized eigenvalue equation (8). Based on the eigenvalues and eigenvectors of the outgoing waves propagating away from the scattering region along the terminal, both the surface Green’s function in Eq. (9) and the group velocity matrix Vm in Eq. (10) are obtained for the m-th terminal in a straightforward calculation. Once the Green’s function for the scattering region GS in Eq. (3) is calculated, the transmittance between the l-th phonon mode in the n-th terminal and the k-th phonon mode in the m-th terminal in Eq. (12) can be calculated based on the matrix tmn in Eq. (11).

The feasibility and validity of the proposed method are demonstrated by two numerical examples. The first chain example is employed with a comparison with the WP simulation. The second example is about a graphene nanojunction with three terminals. The examples illustrate that the improved AGF method is helpful to better understand the heat conduction in complicated structures to the level of individual phonon polarizations.

Reference
[1] Chen X B Liu Y Z Duan W H 2018 Small Methods 2 1700343
[2] Mingo N Yang L 2003 Phys. Rev. 68 245406
[3] Datta S 2005 Quantum Transport: Atom to Transistor New York Cambridge University Press 285
[4] Wang J Wang J S 2006 Phys. Rev. 74 054303
[5] Wang J S Wang J Lu J T 2008 Eur. Phys. J. 62 381
[6] Wang J Wang J S 2009 J. Appl. Phys. 105 063509
[7] Huang Z Murthy J Y Fisher T S 2011 J. Heat Trans. T. ASME 133 114502
[8] Ong Z Y Zhang G 2015 Phys. Rev. 91 174302
[9] Latour B Shulumba N Minnich A J 2017 Phys. Rev. 96 104310
[10] Sadasivam S Waghmare U V Fisher T S 2017 Phys. Rev. 96 174302
[11] Ong Z Y 2018 J. Appl. Phys. 124 151101
[12] Yang L Latour B Minnich A J 2018 Phys. Rev. 97 205306
[13] Zhang L F Wang J S Li B W 2009 New J. Phys. 11 113038
[14] Gu Y F Wu X L Wu H Z 2016 Acta Phys. Sin. 65 248104 in Chinese
[15] Sartipi Z Hayati A Vahedi J 2018 J. Chem. Phys. 149 114103
[16] Sadasivam S Che Y H Huang Z Chen L Kumar S Fisher T S 2014 Ann. Rev. Heat Transfer 17 89
[17] Gu Y F Wang J L 2017 Numer. Heat Tr. B-Fund. 72 71
[18] Zhang W Fisher T S Mingo N 2007 Numer. Heat Tr. B-Fund. 51 333
[19] Gu Y F Wu X L Ni X Y 2016 Numer. Heat Tr. B-Fund. 70 200
[20] Ong Z Y 2018 Phys. Rev. 98 195301
[21] Chen G 2005 Nanoscale Energy Transfer and Conversion New York Oxford University Press 44
[22] Kittel C 2005 Introduction to Solid State Physics New York John Wiley & Sons, Inc 98
[23] Deymier P A 2013 Acoustic Metamaterials and Phononic Crystals New York Springer 16
[24] Gu Y F 2015 Comp. Mater. Sci. 110 345
[25] Saito R Dresselhaus G Dresselhaus M S 1998 Physical Properties of Carbon Nanotubes London Imperial College Press 166
[26] Scuracchio P Costamagna S Peeters F M Dobry A 2014 Phys. Rev. 90 035429